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ABSTRACT 


Context. Recent statistical studies prove that the percentage of RR Lyrae pulsators that are located in binaries or 
multiple stellar systems is considerably lower than might be expected. This can be better understood from an in-depth 
analysis of individual candidates. We investigate in detail the light time effect of the most probable binary candidate 
TU UMa. This is complicated because the pulsation period shows secular variation. 

Aims. We model possible light time effect of TU UMa using a new code applied on previously available and newly 
determined maxima timings to confirm binarity and refine parameters of the orbit of the RRab component in the 
binary system. The binary hypothesis is also tested using radial velocity measurements. 

Methods. We used new approach to determine brightness maxima timings based on template fitting. This can also be 
used on sparse or scattered data. This approach was successfully applied on measurements from different sources. To 
determine the orbital parameters of the double star TU UMa, we developed a new code to analyse light time effect 
that also includes secular variation in the pulsation period. Its usability was successfully tested on CL Aur, an eclipsing 
binary with mass-transfer in a triple system that shows similar changes in the O-C diagram. Since orbital motion would 
cause systematic shifts in mean radial velocities (dominated by pulsations), we computed and compared our model with 
centre-of-mass velocities. They were determined using high-quality templates of radial velocity curves of RRab stars. 
Results. Maxima timings adopted from the GEOS database (168) together with those newly determined from sky surveys 
and new measurements (85) were used to construct an O-C diagram spanning almost five proposed orbital cycles. This 
data set is three times larger than data sets used by previous authors. Modelling of the O-C dependence resulted in 
23.3-year orbital period, which translates into a minimum mass of the second component of about 0.33 OTq. Secular 
changes in the pulsation period of TU UMa over the whole O-C diagram were satisfactorily approximated by a parabolic 
trend with a rate of — 2.2msyr~ 1 . To confirm binarity, we used radial velocity measurements from nine independent 
sources. Although our results are convincing, additional long-term monitoring is necessary to unambiguously confirm 
the binarity of TU UMa. 

Key words, stars: variables: RR Lyrae - binaries: general - methods: data analysis - techniques: photometric - tech¬ 
niques: radial velocities - stars: individual: TU UMa 


1. Introduction 

A significant part of stars are located in double or multiple 
stellar systems. However, reviews of pulsating stars bound 
in binaries (e.g. Szatmary 1990; Zhou 2010) clearly show 
the lack of stellar pairs with an RR Lyrae component. The 
number of currently confirmed binaries comprising an RR 
Lyrae type pulsator can be counted on one hand. 

The binarity of an object can be revealed in many dif¬ 
ferent ways. For example, detection of eclipses, periodic ra¬ 
dial velocity (RV) changes, or regular astrometric shifts in 
a visual binary can serve as a direct proof of binarity. A 
companion of a periodic variable star can also be detected 
indirectly through changes in timings of light extrema, the 
so-called light time effect (hereafter LiTE). RR Lyrae stars 
are generally located at larger distance from Earth, hence 
astrometric detection of binarity is highly unlikely. Since 
the spectra of RR Lyrae stars are influenced by pulsations, 
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discovering the binary nature of stars through changes in 
the position of spectral lines is also difficult (e.g. Fernley & 
Barnes 1997; Solano et al. 1997). Thus the most promising 
methods are detection of eclipses and LiTE. 

In the Large Magellanic Cloud (LMC) three candidates 
for RR Lyraes in eclipsing binaries were detected (Soszyhski 
et al. 2003). However, these objects were identified as op¬ 
tical blends consisting of two objects, RR Lyrae star and 
eclipsing system (Soszyhski et al. 2003; Prsa et al. 2008). 
A very interesting object was identified by Soszyhski et 
al. (2011) and subsequently studied by Pietrzynski et al. 
(2012) and Smolec et al. (2013). This peculiar eclipsing 
system with an orbital period of 15.24 d contains a com¬ 
ponent that mimics an RR Lyrae pulsator. The detailed 
study showed that this object, OGLE-BLG-RRLYR-02792, 
has a very low mass of only 0.26 9TIq ) which is too low for 
classical RR Lyrae stars. Other physical characteristics also 
indicate that this binary component is a very special ob¬ 
ject as a result of an evolution in a close binary followed 
by a life similar to a classical RR Lyrae star. Pietrzynski 
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et al. (2012) included the object in a new class of pulsat¬ 
ing stars called binary evolution pulsators (BEP). OGLE- 
LMC-RRLYR-03541 is another candidate for RR Lyr in 
eclipsing binary (Soszyriski et al. 2009). Nonetheless, the 
similarity in the shape of the light curve and the short or¬ 
bital period of 16.229 d might mean that it belongs to the 
class of BEP (Hajdu 2015, priv. comm.). In addition, it is 
not yet excluded that it could be actually a blend of two 
close stars. 

Szatmary (1990) published a list of various types of pul¬ 
sating stars bound in binary systems that were visually 
identified on the basis of the LiTE in their O-C diagrams. 
However, detailed information and references are missing. 
The catalogue of various binary systems with pulsating 
components from Zhou (2010, version December 2014) con¬ 
tains in the RR Lyrae class TU UMa, OGLE-BLG-RRLYR- 
02792, and several tens of candidates without any closer 
information. Therefore, including some of these objects as 
binary systems (in both catalogues) is at least doubtful. 

Several examples of RR Lyrae stars in binaries were very 
recently identified through the analysis of their O-C dia¬ 
grams. Li & Qian (2014) found that FN Lyr and V894 Cyg 
are probably in pairs with brown dwarves. Hajdu et al. 
(2015) found 20 additional candidates among OGLE bulge 
RRab variables. Liska et al. (2015) presented an analysis 
of 11 new binary candidates located in the Galactic field. 
Nevertheless, all these candidates need to be confirmed 
spectroscopically or using another independent method. 

In this paper, we present a new analysis of the LiTE in 
TU UMa (23-year long cycle, detected by Szeidl et al. 1986) 
that is based on a much wider sample of O-C than was used 
for the LiTE in TU UMa before (Wade et al. 1999). Highly 
accurate photometric observations that cover two thirds of 
the proposed orbital period, are newly available. In Sect. 2 
we briefly discuss the history of observation of TU UMa 
with an emphasis on its binary nature. In Sect. 3 we sum¬ 
marise the characteristics of our data sample. Except for 
data from various sources, we used the original measure¬ 
ments gathered in 2013-2014. We apply our LiTE proce¬ 
dure (described in Sect. 4 and Appendix A) and model the 
TU UMa data in Sect. 5. Other proofs for binarity are dis¬ 
cussed in Sect. 6, and all results are summarised in Sect. 7. 


2. History of TU UMa observations 

TU Ursae Majoris = AN 1.1929 = BD+30 2162 = HIP 
56088 (a = ll h 29 m 48!49, 6 = +30°04'02"4, J2000.0) is a 
pulsating RR Lyrae star of Bailey ab type. According to the 
Variable Star Index 1 (Watson et al. 2006), its brightness 
in V band varies in the range of 9.26-10.24mag (spectral 
type A8-F8) with a period of about 0.558 d. No signs of the 
Blazko effect have been reported for TU UMa up to now. 

The variability of TU UMa was discovered by Guthnick 
& Prager (1929) on Babelsberg plates. Many authors there¬ 
after studied the star using photoelectric photometry and 
spectroscopy. A detailed description of the history of TU 
UMa research was performed by Szeidl et al. (1986). Only 
the most important information about the LiTE is briefly 
mentioned below, since about 180 articles with the keyword 
TU UMa are currently retrievable at the NASA ADS por¬ 
tal. 


1 http://www.aavso.org/vsx/ 


Payne-Gaposchkin (1939) was the first who noted cyclic 
variations in maxima timings and proposed a 12400-day 
(34-year) long cycle. Important results were obtained by 
Szeidl et al. (1986), who mentioned a secular period de¬ 
crease that causes the parabolic trend in the O-C diagram, 
and a probable 23-year (8400 d) variation that is possibly 
caused by the binarity of the star. Saha & White (1990a) 
detected systematic shifts in RVs that indicate binarity, 
but they used very few RV measurements. They modelled 
the LiTE for the first time and determined orbital period of 
7374.5d (20.19yr). Their analysis showed that the proposed 
stellar pair has an extremely eccentric orbit with e = 0.970. 
They considered only a constant pulsation period in their 
model. The effect of neglecting the secular changes of e, 
a sin*, and sin 3 i was tested by Wade et al. (1992). 

The LiTE with secular variation in the pulsation period 
was firstly solved by Kiss et al. (1995), who determined 
more accurate orbital elements and derived an orbital pe¬ 
riod of 8800d (24.1 yr). Wade et al. (1999) collected all 
available maxima timings and also RVs and successfully 
verified results from Saha & White (1990a). They obtained 
five different groups of models of the LiTE with respect 
to different subsets of maxima timings (all maxima with¬ 
out visual values, only photoelectric and CCD values, etc.). 
They derived orbital periods in the range from 20.26 yr 
to 24.13 yr depending on the particular data set and the 
number of fitting parameters (with or without parabolic 
trend). Subsequently, they used nine derived sets of orbital 
elements to reconstruct the orbital RV curve of the pulsat¬ 
ing component and to compare it with shifts in observed 
RVs. Since then, TU UMa has not been studied with regard 
to its LiTE for about 15 years. However, improvements of 
quadratic ephcmeris were performed by Arellano Ferro et 
al. (2013), for instance. Many high-accurate maxima tim¬ 
ings (mainly CCD measurements) were published during 
this 15-year interval. The currently available data exceed 
five proposed orbital cycles. 

3. Data sources 

3.1. GEOS RR Lyrae database 

Since the GEOS RR Lyrae database 2 (Croupe Europeen 
d’Observations Stellaires , Boninsegna et al. 2002; Le 
Borgne et al. 2007) is the most extended archive of times 
of RR Lyrae star maxima, we used this as the main data 
source for our analysis. The version of the database of 
November 2014 contains corrected values for three maxima 
timings from Kiss et al. (1995) that originally contained 
incorrect heliocentric corrections (Wade et al. 1999). 

We discarded all visually determined maxima timings 
because they were inaccurate and omitted two photo¬ 
graphic maxima from Robinson & Shapley (1940) and 
Payne-Gaposchkin (1954), which were replaced by our max¬ 
ima from the DASCH project (see Sect. 3.3). We used also 
the latest maxima timings from Hiibscher (2014); Hiibscher 
& Lehmann (2015) that were not included in the version of 
the GEOS list we used. 

We paid special attention to maxima timings based on 
data from sky surveys such as Hipparcos (ESA 1997) or 
ROTSE = NSVS (Wozniak et al. 2004) for the GEOS val¬ 
ues 3 . These timings were acquired with a special method 

2 http://rr-lyr.irap.omp.eu/dbrr/ 

3 GEOS data marked ‘pr. com.’ were not used. 
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Phase 


Fig. 1 . Differential magnitudes of TU UMa obtained with the 
1-inch telescope folded with the pulsation period (black dots) 
plotted together with the model of V-band observations from 
(Liakos & Niarchos 2011a). 


because in most cases they were determined statistically 
based on the combination of many points that are often 
spread out over a few years. O-C values determined from 
such data sets very often did not follow the general trend of 
O-C dependence 4 . They were therefore omitted. However, 
we re-analysed the original data of these surveys (Sect. 3.3). 
All used maxima timings are available at the CDS portal. 

3.2. Our observations 

As an extension of the GEOS data we also used ten new 
maxima timings gathered by J. Liska in 19 nights between 
December 2013 and June 2014. CCD photometric measure¬ 
ments were performed using two telescopes - three nights 
with a 24-inch Newtonian telescope (vby Stromgren fil¬ 
ters) at Masaryk University Observatory in Brno, and 16 
nights with a small 1-inch refractor (green filter with similar 
throughput as the Johnson V filter, Liska & Liskova 2014) 
at a private observatory in Brno. For the small-aperture 
telescope, five frames with exposure times of 30 s each were 
combined to a single image to achieve a better signal- 
to-noise ratio. The time resolution of such a combined 
frame is about 170 s. The comparison star BD+30 2165 
was the same for both instruments, but the control stars 
were BD+30 2164 (for the 24-inch telescope) and HD 99593 
(for the 1-inch telescope). Maxima timings were determined 
via polynomial fitting 5 or the template fitting described in 
Sect. 3.3. Monitoring with the small telescope resulted in 
the well-covered phase light curve shown in Fig. 1. Except 
for maxima timings determination, the observations were 
also important for detecting possible eclipse (Sect. 6.2). Our 
measurements are available at the CDS portal. 


4 the original value of the maximum timing 2448500.0710 HJD 
from the Hipparcos satellite (Maintz 2005), e.g., has a residual 
value O-C res = 0.0096 d based on model 2 in Sect. 5, but the 
standard deviation of CCD measurements determined from the 
model is only 0.0017 d. 

5 Used for observations with 24-inch telescope where a full 
light curve was not available. 


3.3. Other sources 

We used high-cadence measurements from the Super WASP 
project (Pollacco et al. 2006; Butters et al. 2010) and Pi of 
the Sky project (e.g. Burd et al. 2004; Siudek et al. 2011) 
to determine maximum timings from the individual well- 
covered nights. 

In addition, to extend the O-C dataset as far as pos¬ 
sible, we also analysed data from other large sky surveys 
(Hipparcos, NSVS) and from the project DASCH (photom¬ 
etry from scanned Harvard plates, e.g. Grindlay et al. 2009). 
For these very sparse but very extended data, maxima tim¬ 
ings were estimated using template fitting. The same pro¬ 
cess was applied to the data from the other sources. 

Firstly, we chose the dataset with the best-quality data 
and with the best phase-coverage. Since the data from all 
surveys were of insufficient quality to construct the tem¬ 
plate light curve (e.g. in Hipparcos and NSVS data the 
phase around maximum brightness was not observed), V- 
band measurements from Liakos & Niarchos (2011a) were 
used. We then modelled the shape of the light curve m(t) in 
Matlab with a non-linear least-squares method (hereafter 
LSM) with an n-order harmonic polynomial 

cos ^2 7 rj * p M ° + (fij ^ , (1) 

where Aq is the zeroth level of brightness in mag, Aj are 
amplitudes of the components in mag, t is the time of ob¬ 
servation in Heliocentric Julian Date (HJD), Mq is the zero 
epoch of maximal brightness in HJD, P pu i s is the pulsation 
period in days, and < f>j represents the phase shifts in radi¬ 
ans. After some experiments, we found that a polynomial 
with n = 15 is sufficient for a good template model. 

Nevertheless, using one template curve for various sur¬ 
veys brings some particularities. Firstly, each dataset has 
its own magnitude zero level, and the model has to be scaled 
because different surveys use different filters. Therefore we 
used the whole light curves for our fit. Individual ampli¬ 
tudes Aj and phase shifts <f>j known from the template curve 
remained fixed, but the total amplitude was rescaled using 
the ratio of total amplitudes for the analysed and template 
datasets. The factor was typically close to 1. In addition, 
pulsation period and zero epoch had to be slightly refined 
for each dataset. Subsequently, outliers were iteratively re¬ 
moved. 

Individual observational instruments have different 
spectral sensitivities, therefore it was necessary to verify 
the usability of the same V-band template for different 
datasets. We used photometric measurements in UBVRI 
filters from Liakos & Niarchos (2011a) and compared all 
curves. The main difference is in amplitude, which is highly 
dependent on the mean wavelength (the amplitude de¬ 
creases from 1.3mag in U to 0.6mag in I). When the ampli¬ 
tudes are normalised, the differences between shapes of U, 
B , V, and R curves are almost negligible (comparable with 
the noise level). The same applies for the times of maxima. 
Only /-band light changes differ, apparently. Surveys typi¬ 
cally observed in broad band filters and mean wavelengths 
are not known, or the values were only roughly estimated. 
We can fortunately estimate them using comparison be¬ 
tween amplitudes of survey data and UB VI?/-amplitudes. 
We found that effective wavelengths for all surveys lie be¬ 
tween the wavelengths of the B and R filters; this means 
that our approach (V-band template) is applicable. 


m(t) = A q + ^2 Aj 
j =i 
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Fig. 2. Light curves of TU UMa from different datasets: DA 
- DASCH A, DB - DASCH B, HIPP - Hipparcos, LIAKOS - 
Liakos & Niarchos (20ffa, U-filter), OUR - this paper (1-inch, 
green- filter), WASP - Super WASP (CCD 144), PI - Pi of the 
Sky, NSVS - NSVS, MODEL - used template curves (DASCH 
templates are described in Sect. 3.4). Light curves are vertically 
shifted to better display the light variability. 


Table 1. Numbers of new maxima timings of TU UMa deter¬ 
mined from individual projects and from our observations. 


Hipparcos NSVS 

Pi of the Sky 

DASCH SuperWASP Our 

3 4 

5 

29 21 10 


The good consistence of our template and the observed 
curves was controlled visually (see Fig. 2). For Super WASP 
data, which were obtained using six different CCD cameras, 
a very careful selection of the best nights compared to the 
template was performed to avoid significant trends that are 
present in the data. 

When the data were sparse without well-defined max¬ 
ima, it was necessary to divide the whole dataset into 
smaller subsamples containing typically about 30 points, 
with a time span from several days to hundreds of days. 
The data in these subsamples were then compared with the 
refined template light curve, and the time of maximum was 
determined. With this subsampling we were of course only 
able to estimate the mean time of maximum for the time 
interval of particular subset. However, in comparison with 
the total time span of the O-C values, which cover several 
decades, this method is fully appropriate. The numbers of 
all used maxima timings from particular surveys are listed 
in Table 1. 


In addition, we determined five maxima times from the 
data that were omitted or only poorly determined from the 
analysis in Boenigk (1958), Liakos & Niarchos (2011a,b), 
Liu & Janes (1989), and Preston et al. (1961). 

The uncertainty of each time of maximum determined 
by the template fitting method is influenced mainly by the 
choice of the used model (selection of polynomial degree), 
the time determined for the template maximum, and by 
the quality of the fitted datasets. The first two components 


TU UMa, P = 0.55765 d 



Time [d] 


Fig. 3. Change of the light curve shape caused by long expo¬ 
sure time. The model of the real light curve based on B-band 
observations from Liakos & Niarchos (2011a) (red line) is plot¬ 
ted together with models of observed curves with exposures of 
63min (black crosses) and 200min (blue dots). 


were estimated statistically from the template light curve, 
and their combination was set to 0.0008 d. The influence of 
the third source of uncertainties was individually estimated 
for each light curve directly from the LSM. 

3.4. Photographic measurements - DASCH maxima timings 

Photographic measurements that were obtained with expo¬ 
sure times of an hour and longer are special cases. Light 
curves resulting from these observations differ from real 
curves in amplitude and also in time of minima and max¬ 
ima (for example the models of B-band measurements of 
TU UMa from Liakos & Niarchos (2011a) in Fig. 3). We 
performed a detailed analysis of this problem and will pub¬ 
lish the results elsewhere. Therefore we discuss only the 
most important findings. The extrema timings and ampli¬ 
tude change almost linearly with the duration of the expo¬ 
sure time. Maxima are delayed, minima occurred sooner, 
and the amplitude is lower than for shorter exposures. 

The photographic plates that were digitalised in the 
DASCH project have various integration times, mostly from 
1 to 200 min (Fig. 4), therefore it is very difficult to correctly 
analyse the data (different models of the observing curve 
have to be used). We solved this problem by selecting the 
measurements with similar exposures (they can be approx¬ 
imated with the same template curve; the time difference 
is no more than 0.0005 d). 

At first, we selected the measurements according to the 
exposure - group A (58-68 min, older measurements) and 
group B (30-40 min, newer measurements). Then we calcu¬ 
lated template curves for the mean exposures 63 min and 
35 min (deformed by the exposure) based on measurements 
of Liakos & Niarchos (2011a). These two templates were 
used for the DASCH A and B dataset. 

Our tests showed that for ~ 60 min exposure, for 
instance, the maximum brightness based on exposure- 
deformed light curve is delayed of about +0.004 d in com¬ 
parison with real changes. This difference will be apparent, 
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Fig. 4. Exposure times of the photographic plates of TU UMa 
in the DASCH project. The data were divided into two groups 
according to their integration times - group A with 58-68 min 
(red crosses) and B with 30-40 min (green stars). 


for example, in maxima timings determined using polyno¬ 
mial fitting. Therefore, this shift has to be included in the 
final maxima times. When the template fitting method is 
applied with an improved time-corrected model (template 
from the deformed curve), the exposure-length discrepancy 
is fully reduced to zero. 

4. Modelling the LiTE 

4.1. Light time effect 

The LiTE was suggested at the end of the ninteenth century 
in the Algol system by Chandler (1888), while the first de¬ 
tailed theoretical analysis of the problem was performed by 
Woltjer (1922). Irwin (1952a) solved important equations 
for a part of the direct solution of the LiTE and described 
a graphical way to determine the orbital elements. 

Currently, the LiTE is usually solved very accurately 
by applying equations of motion in a two-body system 
(the equations from Irwin (1952a) included a direct solu¬ 
tion of Kepler’s equation) or using numerical calculations 
of a perturbed orbit in a multiple system. Nevertheless, 
the inverse part of the calculations, in which the best so¬ 
lution is determined (minimisation), still remains a prob¬ 
lem to be discussed - various authors use various meth¬ 
ods, for example, damped differential corrections (Pribulla 
et al. 2000), the LSM (Panchatsaram 1981), the simplex 
(Levenberg-Marquart) method (Wade et al. 1999; Lee et 
al. 2010), or a combination of LSM and simplex method 
(Zasche 2008). For example, the simplex method has sev¬ 
eral versions that differ in the setting of the initial parame¬ 
ters, sorting conditions, or in the size of corrections. Thus, 
obtaining the same results, that is, repeat the process, is 
very difficult, even impossible. To avoid this ambiguity, the 
inverse part of our code was constructed on the basis of the 
non-linear LSM described in (e.g., Mikulasek & Zejda 2013) 
applied to modelling the LiTE (for details see Sect. 4.2 
and Appendix A). This way has previously been applied 
to analyse the LiTE in the AR Aur system (Mikulasek et 
al. 2011; Chrastina 2013) and is similar to the one used by 


Van Hamme & Wilson (2007) and Wilson & Van Hamme 
(2014). 

4.2. Fitting procedure 

The code that we used is written in Matlab. It consists 
of several modules: loading measurements and setting ini¬ 
tial input parameters, direct solution of the LiTE including 
an optional parabolic trend, inverse minimisation method, 
and, finally, selecting the best solution and calculating un¬ 
certainties of individual parameters through bootstrap re¬ 
sampling (Sect. 4.3). 

A prediction of maxima timings T ca i calculated accord¬ 
ing to the relation 


Teal — M 0 + Ppuls x N + A (2) 

contains linear ephemeris, as well as correction A for the 
LiTE. The parameter Mq is the zero epoch of pulsations 
in HJD, Pp U i s is the pulsation period in days, and N is the 
number of pulsation cycle from Mq. 

The correction A for the LiTE includes calculating the 
orbit of a pulsating star around the centre of mass of a bi¬ 
nary in relative units and can be expressed by the equation 
adopted from Irwin (1952a) 


A = A 


(1 c 2 ) sin (V + U) 

7 1 + e cos v 


+ e sin uj 


( 3 ) 


where e is the numerical eccentricity, u> is the argument 
of periastron (usually in degrees), A = a\ sini/173.145 is 
the projection of the semi-major axis of primary (pulsat¬ 
ing) component a± in light days 6 according to the inclina¬ 
tion of the orbit i. and v is the true anomaly. The eccen¬ 
tric anomaly E, which is necessary to determine the true 
anomaly z/, is solved with Kepler’s equation by iteratively 
using Newton’s method with a given precision higher than 
1 x 10 -9 arcsec. Kepler’s equation requires a mean anomaly 
M, which is determined from the orbital period P or bit in 
days and the time of periastron passage To in HJD. 

Another more complex model that includes a parabolic 
trend in O-C (e.g. Zhu et al. 2012), uses the modified 
Eq. (2) in the form of 

Teal = Mo + Ppuls X N + — Ppuls Ppuls x N 2 + A, (4) 

where P pu is is the instantaneous pulsation period at the 
moment Mo and the parameter P pu is = dP pu is/dt is the 
rate of changes in the pulsation period in [dd -1 ]. For easy 
comparison with other RR Lyrae stars, we used prescrip¬ 
tions from Le Borgne et al. (2007). Their parameter 03 = 
1/2 Ppuis Ppuis is the rate of period changes per one cycle in 
[dcycle -1 ], and the rate of period changes /3 in [msd -1 ] is 
/3 = 6.31152 x 10 10 <13 /Ppuis or /3 = 0.07305 x 10 io a 3 /P pu is 
in [dMyr -1 ] 7 . The instantaneous pulsation period at arbi¬ 
trary epoch N is P pu \ s {N) = P pu is(M 0 ) + P pu is x N. 

The first step of the non-linear LSM is linearising the 
non-linear model function (Eqs. 2 or 4) by Taylor decompo¬ 
sition of the first order (see Mikulasek et al. 2006; Mikulasek 

6 the semi-ampl itude of the L iTE changes in O-C diagram is 

then Alitb = A \/l — e 2 * cos 2 uj. 

' Le Borgne et al. (2007) probably used a year with 366 d, 

therefore their constant 0.0732 x 10 10 is slightly different. 
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& Graf 2011) 


a 

leal — T ca i(T, bo) + ^2 A bj 
1=1 


dT ca \(T, b) 


db j 


( 5 ) 


where bj are individual free parameters in vector b. Vector 
bo contains initial estimates of parameters, A bj are their 
corrections, g is the number of free parameters (the length 
of matrix b). 

After linearisation, the problem can be solved in the 
same way as in the linear LSM, but with several neces¬ 
sary iterations to obtain a precise solution. Our code gen¬ 
erates the initial parameters quasi-randomly many times 
from large intervals with limits that are defined by user. 
The derivatives are solved analytically. For more details see 
Appendix A. 

The parameter x 2 (bfc) or its normalised value X^(bfc) = 
X 2 (bfc)/(n — g), where n is a number of measurements, was 
used as an indicator of the quality of the fc-fit. 

Since many maxima timings from the GEOS database 
are given without errors or are often questionable, we used 
an alternative approach to determine the weights. The 
dataset was divided into several groups according to the 
type of observations (photographic, photoelectric, CCD, 
DSLR), and weights were assigned to each of the groups 
with respect to the dispersion of points around the model. 
The weights were improved iteratively. Groups with fewer 
than five points (DSLR) were merged with another group 
with similar data quality to avoid unrealistic weight assign¬ 
ment (CCD+DSLR). During the fitting process, outliers 
differing by more than 5 a from the model were rejected. 
All steps of the analysis were supervised visually. 

The LiTE fitting process does not allow estimating the 
masses of the two stars, but only a mass function /(9Jt) 


_ (®t 2 sin*) 3 _ 4tt 2 (ai sin*) 3 

n + g p 2 bit 


( 6 ) 


where 3Jli, m 2 are masses of the components, i is the in¬ 
clination angle of the orbit, P or bit is the orbital period, and 
ai is the semi-major axis of the primary component. Based 
on the studies of Fernley (1993) and Skarka (2014), we 
adopted as the value for the mass of the RR Lyrae compo¬ 
nent SEJli = 0.55 911© and set the inclination angle to * = 90 0 
(sin* = 1). This allows computing the lowest mass of the 
second component by solving the cubic equation 

m 3 2 - /(an) ml -2 f(m ) soli m 2 - f(m ) m{ = o. (7) 


We expect that the variation in the O-C diagram of 
TU UMa can be well described using Eq. 4 and that other 
possible secular variation of the pulsation period can be ne¬ 
glected (it has a low amplitude or appears to be on a longer 
timescale). 


4.3. Bootstrap-resampling 

The (non-linear) LSM itself gives an estimate of the un¬ 
certainty of all fitted parameters, but the method is ex¬ 
tremely sensitive to data characteristics. A slight change of 
the dataset by adding one single measurement can cause a 
significant difference in the new parameters from the pre¬ 
vious solution. We therefore decided to use a statistical ap¬ 
proach represented by bootstrap-resampling to estimate the 
errors. 


The parameters from the best solution were used as ini¬ 
tial values to fit a new dataset, whose points were ran¬ 
domly selected from the original dataset. This procedure 
was repeated 5000 times. From the scattering of individual 
parameters of these five thousand solutions, we estimated 
their uncertainties. The errors in Tables 2 and 3 correspond 
to 1 er. 


4.4. Test object CL Aurigae: an eclipsing binary with 
probable LiTE and mass transfer 

The code was, among others, tested on the well-known de¬ 
tached binary system CL Aur. This eclipsing binary was 
chosen because it shows the LiTE and a secular period 
change. In addition, CL Aur was studied in a similar way 
three times during the past 15 years (Wolf et al. 1999, 2007; 
Lee et al. 2010). 

The times of minima of CL Aur taken from the O- 
C gateway database 8 (Paschke & Brat 2006) were used 
to construct O-C diagram and to determine parameters 
through the methods described above 9 . Our best model 
(Fig. 5) describes the O-C variations very well in the most 
recent part (precise CCD observations). The old part of the 
O-C diagram with visual and photographic measurements 
is highly scattered, but these measurements were also taken 
into account during the fitting process by assigning them a 
lower weight (model weights for different observation meth¬ 
ods were found in the ratio pg:vis:ccd 1:11.5:403) 10 . We 
conclude that our results are comparable with previous re¬ 
sults (Table 2). This example, as well as additional test¬ 
ing, showed that our code works well and is suitable for 
analysing RR Lyraes with suspected LiTE. The code was 
successfully used to model the LiTE in V2294 Cyg (Liska 
2014). 


5. LiTE in TU UMa and analysis of the 
O-C diagram 

Cyclic changes in O-C diagram of TU UMa have been well 
known for a long time and were analysed in detail by Saha & 
White (1990a) 11 , Kiss et al. (1995), and Wade et al. (1999). 
The parameters determined by these authors are given in 
Table 3. 

Our dataset described in Sect. 3 is much denser and 
more extended than in previous studies (in contrast to 
Wade et al. (1999), who used only 83, we used 253 maxima 
timings). Our O-C values span 113 years, which is mainly 
due to measurements recorded on the Harvard plates (pro¬ 
vided by the project DASCH) in the first half of the twen¬ 
tieth century. Because the complete dataset is not homoge¬ 
neous (it has incomparably better coverage over the last two 
decades with CCD measurements), we analysed LiTE in 
TU UMa in two ways. Model 1 is based on the whole 

8 http://var.astro.cz/ocgate/, all used minima timings are 
available at the CDS portal. 

9 The model was calculated for the half-value of the period. 
Subsequently, results were corrected for this effect. 

10 Wolf et al. (2007) visually distinguished the data quality by 
weights in each category 0, 1, 2 for pg, 0, 1, 2 for visual, and 5, 
10, 20 for CCD observations, respectively. 

11 Several parameters from Saha & White (1990a) e.g. A, a sin i 
were corrected in their erratum (Saha & White 1990b). 
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HJD-2400000 [d] HJD-2400000 [d] 


Fig. 5. O-C diagram of the testing eclipsing binary CL Aurigae 
(double star without an RR Lyrae component) with the 
LiTE and parabolic trend (black circles). The model of changes 
(red line) is based on our parameters from Table 2. 


Table 2. Our determined parameters from the testing object 
CL Aur (right) together with results from previous studies (left). 


Study 

Wolf et al. (2007) 

Lee et al. (2010) 

Our model 

P(Mo) [d] 

1.24437505(18) 

1.24437498(17) 

1.244374881^2 

M 0 [HJD] 

2450097.2712(5) 

2450097.27082(46) 

2450097.27161® 

10 10 P 
[dd- 1 ] 

4.05(6)* 

3.92(55)* 

3-761“ 

10 -10 a 3 
[d cycle -1 ] 

2.52(4) 

2.44(34) 

1° 

CO 

1 + 

4^ -+ 

0 [msyr- 1 ] 

12.8(2)* 

12.4(1.7)* 

11.91® 

0 [dMyr- 1 ] 

0.148(2) 

0.143(20) 

0.137lg° 

Pi [yr] 

21.7(2) 

21.63(14) 

21.6lljg 

To [HJD] 

2443880(80) 

2444072(56) 

24440201^ 

e 

0.32(2) 

0.337(53) 

0.27lg 

u [°] 

209.2(1.2) 

218.9(2.7) 

218lg 

A [light day] 

0.0144(12)* 

0.01378(72)* 

0.01388l2| 

a \2 sin i [au] 

2.49(22)* 

2.38(12) 

2.40l| 

f(m 3 ) [sm©] 

0.034 

0.0290(15) 

0.02971 1 ® 

I< 12 [kms —1 

- 

- 

3.44lg° 


- 

- 

1.04(10) 

Vnin 

144 

198 

203 


Notes. M Parameter calculated using values from the original study. 


dataset 12 , while model 2 describes only photoelectric, CCD, 
and DSLR observations 13 . Since TU UMa experiences sec¬ 
ular period decrease (in Fig. 6 represented by the parabolic 
dot-dashed curve) with the rate l/2P pu i s P pu i s of about 
—2.9 x 1CT 11 days per cycle (Wade et al. 1999), we used 
the complex form of the model (Eq. 4). 

12 Weights for model 1 were found to be of the ratio 
pg : pe: CCD+DSLR 1.0:7.6:14.8, uncertainties are 0.0066 d, 
0.0024 d, and 0.0017 d. 

13 Weights for model 2 were of the ratio pe: CCD+DSLR 
1.00 : 1.38, uncertainties are 0.0020 d and 0.0017 d. 



30000 35000 40000 45000 50000 55000 60000 
HJD-2400000 [d] 

Fig. 6. O-C diagram of TU UMa. Black circles and blue stars 
display the maxima adopted from the GEOS database and new 
maxima determined in this work. The period decrease mani¬ 
fested by the parabolic trend (dot-dashed line) is obvious. Cyclic 
changes due to an orbital motion are also clearly visible. Our 
model of LiTE is represented by the solid red line. The top 
panel shows model 1 with all available data, while the plot in 
the bottom panel shows the situation with only photoelectric, 
CCD, and DSLR observations (model 2). 


Logically, model 1, which is based on the whole data 
set (covering five orbital periods) gives more precise, but 
slightly different results than model 2, which spans only two 
of the 23-year orbital cycles. Nevertheless, high-accurate 
photoelectric and CCD measurements cover the whole or¬ 
bital cycle very well (see Fig. 7). Because of the shorter time 
base, our second model has a lower value of the secular pe¬ 
riod evolution, for example. This happens at the expense of 
an increasing eccentricity (0.66 and 0.69 for models 1 and 
2 , respectively) and change in the argument of periastron 
(181 and 184°). 

In comparison with orbital elements from previous stud¬ 
ies given in Table 3, our results have a higher confidence 
level that is due to the larger and better dataset. Our values 
differ mainly in eccentricity and distance between compo¬ 
nents and in the mass function. All these values were found 
to be significantly lower than those from Saha & White 
(1990a), Kiss et al. (1995), or Wade et al. (1999). Values 
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Table 3. Parameters determined previously (left) and our results (right) for the system TU UMa. The mass limit of the second 
body was estimated from the mass function /(DJI), the orbit inclination (i = 90°), and the mass of the RR Lyrae star DJli = 0.55 DJI© 
adopted from Fernley (1993) and Skarka (2014). Our parameters (right part of the table) were calculated from all maxima timings 
for model 1 and only from photoelectric, CCD, and DSLR measurements for model 2. 


Study 

Saha &; White (1990a) 

Kiss et al. (1995) 

Wade et al. (1999)° 

model 1 

model 2 

PpuiAMo) [d] 

0.5576581097 

0.5576581097 

0.55765817(29)° 

0.557657598120 

0.5576574771!! 

M 0 [HJD] 

2425760.4364 

2425760.4364 

2425760.464(5)° 

2442831.4869lg 

2442831.48641! 

-^puls 

10 -11 [dd -1 ] 

X 

-31.48* 

-10.4(7)* 

—6.99±®g 

—4.551!! 

a 3 = 1/2 -fpuls^puls 
10“ 11 [d cycle -1 ] 

X 

-8.78* 

-2.9(2) 

-1.951® 

-1.271!! 

P = Fpuls [msyr -1 ] 

X 

-9.934* 

-3.3(2)* 

-2.2ll® 

—1.441!! 

P = Rpuls [dMyr- 1 ] 

X 

-0.11498* 

-0.038(3)* 

-0.0255lg° 

—0.01661!! 

Forbit [yd 

20.19* 

24.1(3)* 

23.27(24)* 

23.301® 

23.271! 

T 0 [HJD] 

2425000* 

2447200(50) 

2421585(207) 

2447092 I 4 Q 

24471241!! 

e 

0.970 

0.90(5) 

0.74(10) 

0.6631!° 

0.6861!! 

" H 

196.1* 

178(3) 

183(5)* 

181-311.! 

184.11!;° 

A [light day] 

0.056 

0.023(5)* 

0.0203(35)* 

0.01681! 

0.01721! 

a,i sin i [au] 

10 

4.0(7)* 

3.52(61) 

2.9ll? 0 

2.99!“ 

f(m) [on©] 

- 

0.11(1) 

0.080 

0.046l| 

0.0491! 

®t 2 , min * pn©] 

- 

0.17 

- 

0.3271!! 

0.339!!! 

A'i [kms -1 ] 

60.7* 

11.4(5) 

6.6 

4.971®! 

5.251!° 

x| 

- 

- 

- 

1.05(9) 

1.03(10) 

N max 

~43 

~ 42 

67 

253 

220 


Notes. (*) Parameter calculated using values from the original study, U'J their approach C was selected, 1°) pulsation elements are known 
only from their approach D. 
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Fig. 7. O-C diagram of TU UMa constructed only from pho¬ 
toelectric, CCD, and DSLR measurements after subtracting the 
parabolic trend and phased with the orbital period based on 
model 2. 


from Saha & White (1990a) differ more because they ne¬ 
glected the period decrease. 

The eccentricity is not as extreme as proposed by Saha 
& White (1990a). Its value is comparable with other sys¬ 
tems from Hajdu et al. (2015); Li & Qian (2014). 

Based on our results, it seems that TU UMa is very 
likely a member of a well-detached system with a dwarf 
component with a minimal mass of only 0.33 9 JTq. Since 
no signs of the companion are observed in the light of TU 
UMa, it is probably a late-type main-sequence dwarf star. 


However, we cannot exclude the possibility that it is a white 
dwarf or a neutron star, since we do not know the inclina¬ 
tion. 

Except for the LiTE, which is the most remarkable 
feature of the O-C diagram, changes represented by the 
parabolic trend are also apparent. The progression of the 
dependence suggests secular shortening of the pulsation pe¬ 
riod of TU UMa, which is almost certainly an evolution¬ 
ary effect because the mass transfer, which is responsible 
for period changes in close binaries, can be excluded be¬ 
cause of the very wide orbit of TU UMa. In addition, value 
P = Upuis ~ -2.2 ms yr -1 = -0.026 d Myr” 1 can correspond 
to a blueward evolution of the RR Lyrae component, but 
Le Borgne et al. (2007), for instance, reported a higher me¬ 
dian value /3 = —0.20dMyr _1 for their sample of 21 stars 
with a significant period decrease. 

After subtracting our model 1 from the whole 
dataset (Fig. 8), several photographic measurements (in 
the range of between JD 2425000 and 2429000) are deviate 
more than other values, with a systematic shift of about 
15 minutes (0.01 d). This might indicate that TU UMa un¬ 
dergoes more complex period changes than the LiTE and 
parabolic trend alone. Some possible explanations such as 
cubic trend or an additional LiTE could describe this vari¬ 
ation. However, even though the influence of longer expo¬ 
sure of photographic measurements was taken into account, 
some other instrumental artefact might play a role. 


6. Other proofs for binarity 

Since the LiTE is only an indirect manifestation of binarity, 
it is necessary to prove it in a different way. The analysis 
of the mean RVs can be considered as the most valuable 
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Fig. 8. Residual O-C diagram of TU UMa after subtracting the 
first LiTE model (top panel) and second model - only photoelec¬ 
tric, CCD, and DSLR observations (bottom). Black circles and 
blue stars display maxima adopted from the GEOS database 
and new maxima determined in this work. The jump in the gen¬ 
eral trend of O-C in the range from JD 2425000 to 2429000 (top 
panel) might be an indication of a more complex period change. 

test. In addition to this method, other possible approaches 
with which the binarity of TU UMa can be confirmed are 
discussed in the next sections. 

6.1. Radial velocities 

The known orbital parameters from the analysis of the 
LiTE allow us to predict, but also reconstruct, the RV curve 
from the past. The binarity can then be proved by compar¬ 
ing the model for the orbital RV curve and the spectroscopi¬ 
cally determined centre-of-mass RV. This analysis was first 
performed for TU UMa by Saha & White (1990a) and a 
few years later by Wade et al. (1999). They noted system¬ 
atic shifts in the RV determined at different times. Their 
predictions correlate fairly well with measurements. 

We scanned the literature for RV measurements and 
found nine sources (Table 4). Unfortunately, the last avail¬ 
able dataset with RVs was published in 1997. Saha & White 
(1990a) used only three RV sources, Wade et al. (1999) 
did not give their values. In addition, both authors ignored 


Table 4. Sources of RV measurements for TU UMa, S is a 
source number, AW is a number or RV measurements. 


s 

Publication 

AW 

Lines 

1 

Abt (1970) 

1 

Unknown 

2 

Barnes et al. (1988) 

74 

Metallic 

3 

Fernley & Barnes (1997) 

3? 

OI triplet, (H q ) 

4 

Layden (1993, 1994) 

5 

Hydrogen, Call K 

5 

Liu & Janes (1989) 

60 

Metallic 

6a 

Preston et al. (1961) 

4 

Metallic 

6b 

21 

H 7 , H a , H 8 _n, Ca II K 

7a 

Preston & Paczynski (1964) 

12 

Metallic 

7b 

8+7 

Hydrogen 

8 

Saha Sz White (1990a) 

32 

Metallic 

9 

Solano et al. (1997) 

3? 

Metallic, (H 7 ) 


the RV measurements from Preston & Paczynski (1964). 
Other authors have reported slightly different values de¬ 
termined from the same dataset (Table 5, Col. 2) without 
listing the mean time of observation. Therefore we decided 
to re-analyse all available RV measurements. 

Determining the centre of mass of the RV for a binary 
with a pulsating star is more complicated than for a bi¬ 
nary with non-variable stars (pulsations are often the dom¬ 
inant source of RV changes). Other inconveniences are con¬ 
nected with an RV based on different types of spectral lines 
(e.g. Balmer or metallic lines). They are formed at different 
depths, and therefore RV curves from different lines have 
different shapes, amplitudes, and zero points (shown e.g. by 
Sanford 1949; Oke et al. 1962). Since available RV measure¬ 
ments are based on various lines, it was necessary to unify 
them. This was done using the highly accurate normalised 
template curves from Sesar (2012). Firstly we modelled 
these template curves with an n-order harmonic polyno¬ 
mial. The observed RV curve for the particular spectral line 
was then compared with the polynomial image of the tem¬ 
plate, and the amplitude and the central value of RV curve 
was simultaneously determined by the LSM for all datasets. 
Before this step, measurements were time-corrected for bi¬ 
nary orbit and period shortening (based on model 1), and 
several of the datasets were divided into smaller groups to 
obtain a time resolution of about one year (see Table 5 with 
the determined mean RVs for given epochs corresponding to 
the mean value of observation time). We did not find orig¬ 
inal RV measurements for two studies (Fernley & Barnes 
1997; Solano et al. 1997) and therefore adopted their mean 
RV values and estimated the mean time of observation from 
the information provided in their papers. Finally, the mean 
centre of mass of the RV values were compared with the 
RV model resulting from our LiTE analysis (Fig. 9). The 
points roughly follow the model RV curve. 


An alternative test for binarity using RV curves can be 
performed by comparing the observed RV curves (the top 
panel of Fig. 10) and those in which the orbital RV curve 
from the model is subtracted (the bottom panel of Fig. 10). 
In this figure RV measurements are phased according to the 
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Table 5. Determined values of centre-of-mass velocities for TU 
UMa based on different RV measurements and templates from 
Sesar (2012). The mean values published in different publica¬ 
tions are present for comparison, S is the source number of the 
original RV measurements from Table 4. 


s 

RVp U b 

T m id 

RV our 

errRVour 


[kms -1 ] 

[HJD] 

[kms -1 

[kms -1 ] 

1 

104(35)“ 

2426076 

104“ 

35“ 

2 

90(2) F , 90(2) So 

2443563 

95 

3 



2443941 

98 

3 



2444218 

90 

1 



2444948 

88 

4 

3 

ioi( 3) f , 101(5) 5 ° 

2449520 

101 F 

3 f 

4 

75(17)“, 75(17) F 

2447975 

75“ 

-j^yZ/d 

5 

84.2-“, 84 s “, 84(1)“, 

2446843 

84 

1 


84(1) F , 84(2) So 

2447130 

85 

3 

6a 


2436979 

93 

1 

6b 

92(1) f , S7 h , 

2436648 

93 

2 


93(3) Sa , 92(1 ) f 

2436979 

94 

4 

7a 


2438039 

94 

2 

7b 


2438039 

94 

3 

8 

77 Sa , 77(2) F , 77(2) So 

2446894 

76 

1 

9 

96(3) So 

2449600 

96 So 

3 So 


Notes. Values adopted from ^Fernley & Barnes (1997), 
^Hemenway (1975), ^“■’Layden (1994), (“I Liu & Janes 
(1990), ^Preston et al. (1961), ^“^Saha & White (1990a), 
^ s °^Solano et al. (1997). 


pulsation period 14 . Apparently, the phased RV curve with 
corrected velocities is significantly less vertically scattered 
than without the correction. The residual scatter in the 
bottom panel results from the different metallic lines that 
the RVs were based on. 

Both tests clearly show that TU UMa is very likely 
bound in a binary system. 

6.2. Eclipses 

A detection of eclipses in the light curve (in the appropri¬ 
ate phase of the orbit) would be a strong proof for binarity 
of TU UMa. Several third components of eclipsing binary 
stars that were known only from the LiTE were confirmed 
by detecting additional eclipses (e.g. in the Kepler project, 
Slawson et al. 2011). The probability of catching an eclipse 
in TU UMa is very low because the expected orbital period 
of the binary system is very long and radius and luminos¬ 
ity of the secondary component are probably much smaller 
and lower than for the pulsating star. In addition, the incli¬ 
nation angle of the orbit is unknown. Our two preliminary 
LiTE models allow us to estimate the time of a possible 
eclipse, where the RR Lyrae component should transit the 
secondary one (January - February 2014 or May - June 
2014), but the difference between the two predictions is 


14 The stitching in phase was possible only by taking the 
LiTE and secular period change into account. Without this, the 
curve was scattered horizontally. 


- model 1 model 2 • obs RV - 91 km s _1 



Orbital phase 

Fig. 9. Models of variations in RV caused by orbit of pulsat¬ 
ing component around mass-centre of the binary system (red 
and blue lines) and center-of-mass velocities determined for each 
dataset of RV measurements using template fitting or adopted 
from literature. The visually estimated correction — 91 kms -1 
for systematic velocity mass-centre of the system from Sun ( 7 - 
velocity) was applied. 


too large. However, we attempted to detect the proposed 
eclipse. 

Observations with the small telescope described in 
Sect. 3.2 were dedicated for this purpose. Unfortunately, 
our measurements were insufficient for a reliable decision 
about eclipses. Weather conditions, limited object visibil¬ 
ity and other influences allowed us to observe in only 19 
nights, which is hardly sufficient considering the imprecise 
eclipse prediction. At least we can conclude that no sign 
of an eclipse with an amplitude higher than 0.07mag was 
detected in our data (see Fig. 11). 

7. Summary and conclusions 

We presented a new analysis of a probable LiTE in 
TU UMa. We used published maxima timings from the 
GEOS database (168 values) and added maxima val¬ 
ues from our photometric observations and from the 
SuperWASP and Pi of the Sky surveys. We applied the 
template fitting method to determine maxima from these 
measurements and also from sparse data from the projects 
Hipparcos, NSVS and DASCH. Altogether, we analysed 253 
maxima timing measurements, which is about three times 
more than were used in the dataset in the last study of 
TU UMa by Wade et al. (1999). This large and well cov¬ 
ered dataset allowed us to determine a quadratic ephemeris 
of the pulsations and orbital elements of the binary sys¬ 
tem with much better accuracy than in previous studies 
(Table 3). All analyses were performed with a new code 
written in Matlab that uses a bootstrap method to esti¬ 
mate the errors. We calculated two models: model 1, which 
describes the whole dataset (without visual maxima tim¬ 
ings), and model 2, which describes only high-accurate pho¬ 
toelectric, CCD and DSLR maxima. The second model is 
based on data with a significantly shorter time span than 
for model 1. 
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Corrected phase 


Fig. 10. Radial velocity curves from the metallic lines of TU 
UMa from different publications phased with the pulsation 
period corrected for the LiTE and secular period changes. 
Uncorrected observed RVs (top) and corrected values after sub¬ 
tracting the changes in RV caused by orbital motion based on 
our model 1 (bottom). RV values corrected for binary orbit are 
evidently less scattered than uncorrected RVs. 


Fig. 11. Residuum of the light curve of TU UMa after subtract¬ 
ing the harmonic polynomial model. No signs of an eclipse with 
an amplitude higher than 0.07 mag was detected in the green 
band. 

eclipse is highly unlikely (wide orbit, unknown inclination, 
and other important parameters) we attempted to detect 
them, but were unfortunately not successful. 

Binarity is expected manifest itself in cyclic changes of 
the mean RV. We adopted RV measurements from nine in¬ 
dependent sources and corrected their values according to 
our model by subtracting the LiTE and secular changes. 
When an observed RV curve was phased with the pulsation 
period, we obtained the typical RV curve for RR Lyrae, 
which was scattered. The scatter significantly decreased 
when our model was applied. 

We also determined central RV values for each RV 
dataset using pulsation templates for different spectral lines 
from Sesar (2012). We compared these values with our 
model of orbital RV variations based on orbital parameters 
known from the LiTE. The correlation is evident (Fig. 9). 

The two successful proofs are important for confirming 
the binarity of TU UMa. However, only long-term spectro¬ 
scopic measurements covering the whole orbital cycle can 
unambiguously confirm that TU UMa is indeed a member 
of a binary system. 


The second model gives a lower value of the period- 
decrease rate (/? = -P pu is ~ —1.4msyr _1 ), which causes 
the eccentricity to become higher (e ~ 0.69) than in the 
first model (/3 ~ — 2.2msyr~ 1 , e ~ 0.66). For comparison, 
Arellano Ferro et al. (2013) give j3 = — 1.3msyr _1 with¬ 
out fitting the LiTE. Nevertheless, both our models have 
a lower eccentricity value, semi-major axis of pulsating 
component a i sin* (2.9au or 3.0au), and semi-amplitude 
of RV variations of the pulsating star K\ (5.0kms _1 or 
5.3kms _1 ) than in previous works. Our values of the or¬ 
bital period (identically 23.3 yr) and argument of periastron 
u> (181° or 184°) are comparable with values determined 
by previous authors. In addition, the lowest mass limit of 
the secondary component (0.33 ©Iq or 0.34 97t 0 ) was deter¬ 
mined with an assumption for the mass of the RR Lyrae 
component of 0.559110. 

The binary nature was tested in several ways. Firstly, 
our models of the orbit gave predictions of possible eclipses. 
Although the prediction was highly inaccurate and an 


Acknowledgements. The DASCH project at Harvard is grateful for 
partial support from NSF grants AST-0407380, AST-0909073, and 
AST-1313370. This paper makes use of data from the DR1 of the 
WASP data (Butters et al. 2010) as provided by the WASP consor¬ 
tium, and the computing and storage facilities at the CERIT Scientific 
Cloud, reg. no. CZ.1.05/3.2.00/08.0144, which is operated by Masaryk 
University, Czech Republic. This research has made use of NASA’s 
Astrophysics Data System. Work on the paper has been supported 
by LH14300. MS acknowledges the support of the postdoctoral fel¬ 
lowship programme of the Hungarian Academy of Sciences at the 
Konkoly Observatory as a host institution. We are very grateful to 
the anonymous referee, who significantly helped to improve this pa¬ 
per. 


References 

Abt, H. A. 1970, ApJS, 19, 387 

Arellano Ferro, A., Pena, J. H., & Figuera Jaimes, R. 2013, Rev. 
Mexicana Astron. Astrofis., 49, 53 

Barnes, T. G. III., Frueh, M. L., Moffett, T. J. et al. 1988, ApJS, 67, 
403 

Boenigk, T. 1958, Acta Astron., 8, 13 


11 































Liska et al.: Light time effect in TU Ursae Majoris 


Boninsegna, R., Vandenbroere, J., Le Borgne, J. F., Sz Geos Team 
2002, IAU Colloq. 185: Radial and Nonradial Pulsations as Probes 
of Stellar Physics, 259, 166 

Burd, A., Cwiok, M., Czyrkowski, H., et al. 2004, Astronomische 
Nachrichten, 325, 674 

Butters O. W., West, R. G., Anderson, D. R. et al. 2010 A&A, 520, 
10 

Chandler, S. C. 1888 AJ, 7, 165 

Chrastina, M. 2013, Study of fast light changes of interacting binaries, 
Ph.D. Thesis, DTPA of Masaryk University, Brno, Czech Republic 
ESA 1997, The Hipparcos and Tycho Catalogs, SP-1200 
Fernley, J. 1993, A&A, 268, 591 
Fernley, J., Sz Barnes, T. G. 1997, A&AS, 125, 313 
Grindlay, J., Tang, S., Simcoe, R. et al. 2009, ASPC, 410, 101 
Guthnick, P., Sz Prager, R., 1929, Beob. Zirk., 11, 32 
Hajdu, G., Catelan, M., Jurcsik, J., et al. 2015, MNRAS, 449, L113 
Hajdu, G., 2015, priv. comm. 

Hemenway, M. K. 1975, AJ, 80, 199 

Hiibscher, J. 2014, Information Bulletin on Variable Stars, 6118, 1 
Hiibscher, J., Sz Lehmann, P. B. 2015, Information Bulletin on 
Variable Stars, 6149, 1 
Irwin, J. B. 1952a, ApJ, 116, 211 
Irwin, J. B. 1952b, ApJ, 116, 218 

Kiss, L. L., Szatmary, K., Gal, J., Sz Kaszas, G. 1995, Information 
Bulletin on Variable Stars, 4205, 1 

Layden, A. C. 1993, The Metallicities and Kinematics of the Local 
RR Lyrae Variables, Ph.D. Thesis, Yale University, New Haven, 
USA 

Layden, A. C. 1994, AJ, 108, 1016 

Le Borgne, J. F., Paschke, A., Vandenbroere, J., et al. 2007, ASzA, 
476, 307 

Lee, J. W., Kim, Ch.-H., Kim, D. H., et al. 2010, AJ, 139, 2669 
Li, L.-J., Qian, S.-B., 2014, MNRAS, 444, 600 

Liakos, A., Sz Niarchos, P., 2011a, Information Bulletin on Variable 
Stars, 6099, 1 

Liakos, A., Sz Niarchos, P., 2011b, Information Bulletin on Variable 
Stars, 5990, 1 

Liska, J., Liskova, Z., 2014, Information Bulletin on Variable Stars, 
6124, 1 

Liska, J., 2014, Information Bulletin on Variable Stars, 6119, 1 
Liska, J., Skarka, M., Zejda, M., Sz Mikulasek, Z. 2015, 
arXiv: 1504.05246 

Liu, T., & Janes, K. A. 1989, ApJS, 69, 593 
Liu, T., Sz Janes, K. A. 1990, ApJ, 354, 273 
Maintz, G. 2005, A&A, 442, 381 

Mikulasek, Z., Wolf, M., Zejda, M., Sz Pecharova, P. 2006, Ap&SS, 
304, 363 

Mikulasek, Z. Sz Graf, T. 2011, arXiv: 1111.1142 
Mikulasek, Z., Ziziiovsky, J., Zejda, M., et al. 2011, in Magnetic 
Stars, Proceedings of the International Conference, held in the 
Special Astrophysical Observatory of the Russian AS, August 27- 
September 1, 2010, Eds: I. I. Romanyuk and D. O. Kudryavtsev, 
431 

Mikulasek, Z., Sz Zejda, M. 2013, in Uvod do studia promennych 
hvezd, ISBN 978-80-210-6241-2, Masaryk University, Brno 
Oke, J. B., Giver, L. P., Searle, L. 1962, ApJ, 136, 393 
Paschke, A., Brat, L. 2006, OEJV, 23, 13 
Panchatsaram, T. 1981, Ap&SS, 77, 179 

Payne-Gaposchkin, C., 1939, Proceedings of the American 

Philosophical Society, 81, 189 

Payne-Gaposchkin, C. 1954, Annals of Harvard College Observatory, 
113, 151 

Pietrzynski, G., Thompson, I. B., Gieren, W., et al. 2012, Nature, 
484, 75 

Pollacco, D. L. , Skillen, I., Collier Cameron, A., et al. 2006, PASP, 
118, 1407 

Preston, G. W., Spinrad, H. Sz Varsavsky, C. M. 1961, ApJ, 133, 484 
Preston, G. W. Sz Paczynski, B. 1964, ApJ, 140, 181 
Pribulla, T., Chochol, D., Tremko, J. et al. 2000, Contributions of 
the Astronomical Observatory Skalnate Pleso, 30, 117 
Prsa, A., Guinan, E. F., Devinney, E. J., Sz Engle, S. G. 2008, A&A, 
489, 1209 

Robinson, L. V., Sz Shapley, H. 1940, Annals of Harvard College 
Observatory, 90, 27 

Saha, A., & White, R. E. 1990a, PASP, 102, 148 
Saha, A., Sz White, R. E. 1990b, PASP, 102, 495, Erratum 
Sanford, R. F. 1949, ApJ, 109, 208 
Sesar, B. 2012, AJ, 144, 114 


Siudek, M., Malek, K., Mankiewicz, L., et al. 2011, Acta Polytechnica, 
51, 68 

Skarka, M. 2014, MNRAS, 445, 1584 

Slawson, R. W., Prsa, A., Welsh, W. F., et al. 2011, AJ, 142, 160 
Smolec, R., Pietrzynski, G., Graczyk, D., et al. 2013, MNRAS, 428, 
3034 

Solano, E., Garrido, R., Fernley, J., Sz Barnes, T. G. 1997, A&AS, 
125, 321 

Soszynski, I., Udalski, A., Szymanski, M., et al. 2003, Acta Astron., 
53, 93 

Soszynski, I., Udalski, A., Szymanski, M. K., et al. 2009, Acta Astron., 
59, 1 

Soszynski, I., Dziembowski, W. A., Udalski, A., et al. 2011, Acta 
Astron., 61, 1 

Szatmary, K. 1990, Journal of the American Association of Variable 
Star Observers, 19, 52 

Szeidl, B., Olah, K., Sz Mizser, A. 1986, Communications of the 
Konkoly Observatory Hungary, 89, 57 
Van Hamme, W., Sz Wilson, R. E. 2007, ApJ, 661, 1129 
Wade, R. A., Saha, A., White, R. E., Sz Fried, R. 1992, Bulletin of 
the American Astronomical Society, 24, 1225 
Wade, R. A., Donley, J., Fried, R., et al. 1999, AJ, 118, 2442 
Watson, C. L., Henden, A. A. Sz Price, A. 2006, Society for 
Astronomical Sciences Annual Symposium, 25, 47 
Wilson, R. E., Sz Van Hamme, W. 2014, ApJ, 780, 151 
Wolf, M., Sarounova, L., Broz, M., Horan, R., 1999, Information 
Bulletin on Variable Stars, 4683, 1 
Wolf, M., Kotkova, L., Brat, L., et al. 2007, Information Bulletin on 
Variable Stars, 5780, 1 

Woltjer, Jr., J. 1922, Bull. Astron. Inst. Netherlands, 1, 93 
Wozniak, P. R., Vestrand, W. T., Akerlof, C. W., et al. 2004, AJ, 
127, 2436 

Zasche, P. 2008, arXiv:0801.4258 
Zhou, A. -Y. 2010, arXiv: 1002.2729v6 

Zhu, L.-Y., Zejda, M., Mikulasek, Z., et al. 2012, AJ, 144, 37 


Appendix A: Application of the non-linear LSM on 
calculating the LiTE 

We assumed a group of n maxima timings given in 
Heliocentric Julian Date 15 (HJD) determined from obser¬ 
vations. Times are inserted in the column vector y with a 
size nxl. Each time of an observed maximum T) has a cor¬ 
responding uncertainty cq. We assumed that the quality of 
the /-measurement can be quantified by the /-weight using 
the relation wi = er ; ~ 2 , and these weights were inserted in 
the column vector w. For a correct calculation, the weights 
were normalised (the average value of weights is w = 1) 
and were inserted in a square matrix W = diag(w). 

In the next step, we selected the equation for the model 
function T ca i(T, b) with unknown parameters in vector b . 
Changes in the position times of the maxima for a pulsating 
star that are caused by the LiTE (periodic changes in the 
O-C diagram) can be expressed by 

T ca i(T,b) =M 0 + 

P puls xN + A, (A.l) 

where Mq is the zero epoch of pulsation in HJD, P pu i s is 
the pulsation period in days, N is the number of pulsation 
cycles from Mo, and parameter A is the correction for the 
LiTE. The integer number of the pulsation cycle (epoch) 
N is 

N = round ■ (A.2) 

15 Times of maxima in Barycentric Julian Date should be used, 
but the barycentric correction is lower than the accuracy of max¬ 
ima times from the GEOS database, which contains times in 
HJD valid to the fourth decimal place. The accuracy of the de¬ 
termined maxima is likewise mostly lower - especially for pho¬ 
tographic or visual measurements. 
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An optional more complex model that includes the and then the sum is 
parabolic trend in the O-C diagram (e.g. Zhu et al. 2012) 
uses the modified Eq. (A.l) in the form of 


x 2 (b) = ^5rLd,i = E 


Teal(T, b) — M 0 + Ppuls X N + — Ppuls Ppuls x N 2 + A, 

(A.3) 


Ti — TAu (Ti , b) 


cn 


= £[«?«*]■ 


(A.12) 


where P pu i s is the instantaneous pulsation period at the 
moment Mo, and parameter P pu i s = dP pu i s /dt is the rate 
of changes in the pulsation period. 

The correction A for the LiTE includes calculating the 
orbit of the pulsating star around the binary mass centre 
in relative units and is given by the equation adopted from 
(Irwin 1952a) 


A = A 


2 sin(u + w) 
1 + e cos v 


+ e sin uj 


(A.4) 


where e is the numerical eccentricity, v the true anomaly, uj 
the argument of periastron in degrees, and A is a constant 
in light days, which compares the shift in a radial position 
to the time delay caused by the constant speed of light. The 
true anomaly v is calculated from 


v 1 + e E 

tan— = \ -tan—, (A.5) 

2 V 1 - e 2 v ' 

and the eccentric anomaly E is iteratively determined in 
our code by Newton’s method from Kepler’s equation 


E = M + e sinE. 


The mean anomaly M is in the form 


M = 


2 tt(T —Tp) 

I orbit 


(A. 6 ) 


(A. 7 ) 


where the orbital period P 0 rbit is in days and the time of 
periastron passage To in HJD. 

The constant A is the projection of the semi-major axis 
of the pulsating component a\ on the unit light day 


di sini au ai sin I 

A = —-= —--, (A. 8 ) 

86400 c 173.145’ v ’ 

where i is the inclination angle of the orbit in degrees, au 
is length of the astronomical unit in metres, c is the speed 
of the light in vacuum in ms^ 1 . The semi-amplitude of 
LiTE changes in the O-C diagram in days is then 


Its normalised value x^(bfc) = y 2 (b k)/{n — g), where n 
is the number of measurements and g is the number of 
free parameters (the length of matrix b), was used as an 
indicator of the quality of the k- fit. 

The first step of the non-linear LSM is linearisation of 
the non-linear model function (Eqs. A.l or A.3) by a Taylor 
decomposition of the first order (see Mikulasek et al. 2006; 
Mikulasek & Graf 2011) 


Teal = T ca l(T, b 0 ) + £} A bj b()) 

1=1 




(A.13) 


where bj are individual free parameters from the vector b, 
vector bo contains the initial estimates of the parameters, 
and Abj are their corrections in the vector Ab. 

Important equations for solving non-linear LSM by ma¬ 
trices are 


U = X T W Ay, V = X T WX, 

H = V -1 = (X T WX) _1 , Ab = HU, 


(A. 14) 


from which the new parameters in vector bb can be calcu¬ 
lated as bi = b 0 + Ab. The difference between observed 
values and model is Ay — y T C ai- The matrix with the 
derivatives has the form 


<9T cal dTcai <9T cal dTAi <9T ca i 
dM 0 ’ <9 -Ppuls ’ <9To ’ dPorbit ’ dA ’ 
dPcal 9T C a\ 

’ duo ’ de 


(A.15) 


where individual derivatives are also presented 

9T ca i <9T cal OA 8u dE dM 

dM 0 ’ dT 0 du dE dM dT 0 ’ 

dT C ai dT ca i dA du dE dM 

dPpuis ’ dPorbit du dE dM dP olh it ’ 

dTcai (1 - e 2 ) sin(i' + cj) 

= -T5- he smu; ’ 

dA 1 + e cos u 


Alite = A x/T^ - e?~cos?Lu. (A.9) 

Subsequently, the observed values of the time of maxi¬ 
mum can be compared with those from the model obtained 
from Eqs. A.l and A.3. Their difference for a given set of 
parameters is equal to 


dJA i 
dui 


= A 


(1 — e 2 ) cos(j/ + w) 
1 + e cos u 


+ e cos oj 


dT ca \ dA dA du dA du dE 

de de ^ du de ^ du dE de ' 


STi = T; — T cali ;(T], b). (A.10) 

The LSM described in detail in (Mikulasek & Zejda 2013) 
points that the best model has the lowest sum of squares 
of residuals between observation and model. The modified 
form of the LSM used the weighted form 


Other necessary derivatives are 


dA _ A (1 — e 2 ) 
du 1 + e cos u 


cos (u + uj) + 


e sin)!' + uj) sin u 


1 + e cos u 


^T m od,l — — 


STi T,-T caM (r,, b) 




<?l 


(AT!) f=A 


sin);' + uj) [2 e + (1 + e 2 ) cos u] 
(1 + e cos u ) 2 
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( rna u \ 

2 

dv 11 + e 

cos- 

du sin v 

dE V 1 -e 

1~|J 

de 1 — e 2 


dE _ 1 dE 

dM 1 — e cos E ’ de 


sin E 

1 — e cos E ’ 


dM _-2tt(T-To) dM _ -2n 

dPorbit ^orbit ’ 9To Porbit 


Matrix X will be expanded by about one additional member 
to calculate the parabolic trend according to Eq. A.3 


<9T ca i (9T ca i dT ca i dT ca \ dT ca \ 

M^’dT^ , ~^’dTu~ t , ~dT’ 

dT caX dT ca \ dT ca i 
’ duj ’ de c)Pp U i a J ’ 


(A.16) 


where two of X members are in the form 


= N + - P pu i s x N 2 

dPp U i s 2 puls 


5Tcal = i Pp U l s x N 2 


dR 


puls 


The determined parameters allow calculating the radial 
velocity (RV) changes caused by the secondary component 
(e.g. Irwin 1952b) 

RV i = 7 + Ki [cos(> + oj) + e cosw], (A.17) 


where 7 is systematic velocity mass-centre of the binary 
system from the Sun in kms ' 1 ( 7 -velocity) and Ki is the 
semi-amplitude of RV changes in kms -1 given by 


K x 


2 7 r ai sin i au 

8.64 X 10 7 Porbit\/(l — e 2 ) ’ 


(A.18) 


where the projection of the semi-major axis a\ sini is in 
au, the constant au in meters, and the orbital period P or bit 
in days. 
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